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We explore a minimal paradigm for thermalization, consisting of two weakly-coupled, low dimen- 
sional, non-integrable subsystems. As demonstrated for Bose-Hubbard trimers, chaotic ergodicity 
results in a diffusive response of each subsystem, insensitive to the details of the drive exerted on 
it by the other. This supports the hypothesis that thermalization can be described by a Fokker 
Plank equation. We also observe, however, that Levy-flight type anomalies may arise in mesoscopic 
systems, due to the wide range of time scales that characterize 'sticky' dynamics. 



The emergence of irreversibility from reversible Hamil- 
tonian mechanics remains an open fundamental question, 
even after a century of effort. Recent advances in com- 
putational as well as experimental technique may at last 
bring answers within reach. The biggest challenge of this 
quest is the sheer technical difficulty of solving the Hamil- 
tonian evolution of quantum many-body systems, even 
when they are quite small and isolated. In this Letter 
we propose to leap over a significant barrier of under- 
standing, by using Hamiltonian results from a tractable 
but non-trivial system, to support an extension of an es- 
tablished phenomenological theory, into a substantially 
more challenging regime. The result we thereby derive is 
a simple theory that can then both guide, and be tested 
by, subsequent numerical investigations, as well as cur- 
rently feasible experiments. 

We address the thermalization of two nonlinear Hamil- 
tonian subsystems that are weakly coupled together, 
where the combined system is isolated and undriven. 
The equilibration of such subsystems is postulated in 
the Zeroth Law of Thermodynamics, reflecting the as- 
sumption that microscopic dynamics is unobservably 
fast, while slower macroscopic dynamics remains non- 
trivial. Accordingly, weakly coupled subsystems, each 
having strong internal interactions, provide the mini- 
mal paradigm for the emergence of thermodynamics from 
closed-system mechanics. 

Following the Fermi-Pasta-Ulam numerical experi- 
ment, most studies of dynamical equilibration have his- 
torically focused on large, extended systems [11121, where 
the treatment of even one strongly interacting system is 
quite impossible in microscopic detail. With experimen- 
tal access to controlled mesoscopic systems, attention has 
more recently been drawn to thermalization phenomena 
in small systems, taking into account dynamical chaos 
[31 Hj and quantum effects [IHI] ■ The traditional anal- 
ysis of thermalization has nonetheless largely remained 
within the assumptions inherited from the macroscopic 
problem. It is common to assume that at least one of 
the two coupled systems is "big" , and hence can be dras- 
tically approximated, either as a phenomenologically de- 
scribed reservoir, or as a time-dependent external param- 



eter. The present Letter is motivated by the realization 
that the study of isolated thermalization of two subsys- 
tems is no longer so unthinkably intractable. It is merely 
extremely difficult. Our proposal is to leverage our un- 
derstanding of driven chaotic systems to overcome this 
difficulty, by viewing each subsystem as driving the other. 

The statistical approach.— The statistical descrip- 
tion of driven chaotic systems by means of a Fokker- 
Planck equation (FPE) for their energy distribution [S]- 
[T3] is based on the ergodic adiabatic theorem [T3] . Quan- 
tum and classical systems can be embraced in a unified 
notation by writing the energy e as a function of the 
phase space volume n of the constant-energy hypersur- 
face. The density of states is g{e) = dn/de, and the 
micro-canonical inverse temperature is /3(e) — d\n{g) /de. 
Upon quantization, the Wigner-Weyl formalism implies 
that n corresponds to the discrete index of the energy 
levels e„, and if these levels are dense enough, they 
can be approximated as a quasi-continuum. One then 
makes a coarse-grained description of the slow evolution 
of the system, and derives an FPE to describe the evo- 
lution of the time-dependent energy probability distribu- 
tion p{£, t). 

FPE for a driven system.— If a chaotic system is 
driven weakly, its energy changes slowly, and p{e, t) obeys 
a probability-conserving FPE, whose diffusion term has 
a coefficient D, proportional to the strength of the driv- 
ing [SHTOl [T5] . By Liouville's theorem, a distribution 
p{e) oc g{e) should be a time-independent solution of the 
FPE. Hence it is deduced that the drift term in the FPE 
is universally related to D, and the complete phenomeno- 
logical equation is established. 

FPE for coupled subsystems. We now extend 
the single-system FPE phenomenology [SUIO] to the 
case of thermalization of two subsystems. Each sub- 
system (i = 1, 2) is characterized by its density of states 
gi{£i), and by its microcanonical inverse temperature Pi. 
Thanks to conservation of energy the thermalization is 
within subspaces of constant energy ei(ni) +£2{n2) — E. 
Accordingly we set ei = e, and £2 = E — e, and construct 
an FPE for the probability density p{£, t) that describes 
how the energy is divided between the two subsystems. 
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It again follows from Liouville's theorem that an er- 
godic distribution p{e) oc g{e) = g\{e)g2{£ — e), should 
be a stationary solution. This fixes the form of the FPE, 
and implies the functional form of the drift term: 
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It is important to notice that the diffusion coefficient D 
may depend on e. The optional way Eq.([2| of writing this 
FPE demonstrates that the 'drift velocity' A is related 
to the diffusion as follows: 
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To see more clearly the connection of Eq. ^ with tradi- 
tional thermodynamics, assume that each of the subsys- 
tems is prepared independently in a canonical state, with 
temperature Ti, such that gi{ei) exp(— e^/Ti) describes its 
energy distribution. Integrating both sides of Eq.(|3| with 
this probability measure, and integrating by parts the 
first term on the right, one obtains [ij a mesoscopic Ein- 
stein relation, like those previously derived [T5l [HI using 
master-equation or fluctuation-theorem approaches: 
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This result offers insight into the distinct behaviors of 
micro-canonical energy fluctuations and canonical aver- 
ages. The canonical version Eq. Q implies that energy 
always flows from the higher to the lower canonical tem- 
perature, but the more general mesoscopic version Eq.([3| 
implies that energy flow is not necessarily from the higher 
to the lower micro-canonical temperature, and may de- 
pend on the functional form of D{e). This is not in con- 
tradiction with the Zeroth Law of thermodynamics: en- 
ergy fluctuates between finite systems in equilibrium, and 
so their micro-canonical temperatures need not be equal. 
The ergodic solution p oc g{e) around which Eq. ([T]) has 
been constructed implies only that the most probable e 
is the one for which /3i(e) = /32{£ — £)■ 

Fluctuation-dissipation phenomenology. The 
derivation of Eq. ([ij is phcnomenological, but based on 
simple assumptions that can be tested. Since these in- 
clude weak coupling, it is further consistent to compute 
D using the Kubo formula. Writing the interaction as 
H — Q'-^'Q''^-', and defining S^'^*'(a;) as the power spec- 
trum of the fluctuating variable it reads [17] 
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With this addition, the FPE phenomenology provides a 
generalized fluctuation-dissipation relation that connects 



the systematic energy flow between the subsystems with 
the intensity of the fluctuations. 

Reasoning.— When two undriven subsystems are 
coupled to each other, the effect of one subsystem (call 
it "agent") on the other (call it "system") is like that 
of driving. For the purpose of obtaining Eq. ([I]) we have 
assumed that the interaction results in diffusion that can 
be calculated using Eq. ([5| . Future studies of coupled 
systems must test this assumption in full, but one key 
point remains to be established with regard to the driven 
single-subsystem dynamics: The agent-system interac- 
tion will typically couple many quantum levels, even if it 
is weak in classical terms; In such strongly non-adiabatic 
circumstances, energy diffusion, and the applicability of 
Eq.([5]), have not been demonstrated Below we com- 
plete this paper by a numerical demonstration that high- 
lights the role of chaos in obtaining diffusive dynamics 
for a driven subsystem, supporting the feasibility of the 
above reasoning for experimentally relevant systems. 

Testing ground.— Few-mode Bose-Hubbard systems 
are a promising testing ground, since they are experimen- 
tally accessible and highly tunable [TH], and also theo- 
retically tractable by a wide range of techniques. Since 
boson number is conserved, their Hilbert spaces are of 
finite dimension, and yet their classical dynamics can be 
non-integrable. The smallest Bose-Hubbard system ad- 
mitting chaos without external driving is the three-mode 
trimer P^H^ , described by the Bose-Hubbard Hamilto- 
nian (BHH): 
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Here i = 0,1,2 label the three modes, and aj are 
canonical destruction and creation operators in second 
quantization, K is the hopping frequency, and U is the 
on-site interaction. The Hamiltonian % commutes with 
the total particle number M — a\ai, and hence, with- 
out loss of generality, we regard N as having a defi- 
nite value N . Driving is then implemented by setting 
K = Kq + KdSm{nt). Consequently the total Hamilto- 
nian has the structure Ho + f{t)W, where the perturba- 
tion operator W is identified as the first sum in Eq.(|6|, 
and the driving field is f{t) = {Kd/2) sm{nt). 

Chaoticity.— The underlying classical dynamics is 
defined [aj by replacing the operators a.; in the Heisen- 
berg equations of motion |26j with complex c-numbers 
riie^'^K In the absence of driving, up to trivial rescal- 



ing, the classical equations depend only on the single di- 
mensionless parameter u = NU/Kq. The chaoticity of 
the motion that is generated by Ho is reflected in the 
local level statistics, and can be quantified by the Brody 
parameter 0<q<l |27], such that q=0 indicates a Pois- 
sonian level-spacing distribution (characteristic of inte- 
grable dynamics), while higher values indicate the ap- 
proach to Wigner level-spacing distribution (indicating 



FIG. 1: The energy spectrum of the unperturbed BHH. In 
panel (a) the scaled eigen-energies £„ of Ho are plotted versus 
the scaled interaction parameter u, for A'^ = 35 particles. The 
level spacing statistics is characterized by the Brody param- 
eter (0 < g < 1), which is displayed in Panel (b) for a system 
with A'^ = 120 particles. In the energy range where the motion 
is chaotic g ~ 1. Square symbols indicate the preparations 
that were used for the simulations in Figj2] 



chaotic dynamics). 

Fig[T^ displays the spectrum e„ , obtained by numerical 
diagonalization of Hq, as a function of u. For graphical 
presentation we shift and scale the energy spectrum, for 
each u, into the same range e G [0, 1], such that e = and 
£ = 1 are the ground energy Eq and the highest energy 
^'max respectively. By plotting q vs {u,e), as in Fig|T|D, 
we can identify the e range within which the motion is 
chaotic at any given value of u. See [ij for technical de- 
tails. We have verified the implied chaoticity by plotting 
representative classical Poincare cross-sections. 

In the numerical simulations we consider an = 50 
particle system with two representative values of u. The 
case u = 5, for which there is a wide chaotic range 
0.2 < e < 0.6, is contrasted with u = 50, for which the 
motion is globally quasi-integrable due to self-trapping. 

Energy diffusion.— In the absence of driving the en- 
ergy is a constant of motion. Driving induces transitions 
between energy eigenstates, leading to a time-dependent 
spread in energy Ae(t). This dispersion is defined as 
the square root of the variance Var(e„) that is associated 
with the probability distribution 



Pn{t) = (£n|*(t)) 
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In Fig[2]we plot the time evolution of the quantum energy 
distribution Pn{t) in response to driving that is quantum 
mechanically large (many levels are mixed) but classi- 
cally small {Kd <C Kq). We contrast the response in the 
chaotic (u = 5) and in the quasi-integrable (u = 50) 
regimes. Dramatic differences are observed. In both 
cases, the energy distribution in the very early stages 
of the evolution reflects the band profile of the pertur- 
bation matrix Wn^noj where rig is the initial level, as ex- 
pected from time-dependent first-order perturbation the- 
ory. Later in the evolution higher orders of perturbation 
theory dominate. This leads in the quasi-integrable case 
to Rabi-like oscillations that have no relation to the clas- 
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FIG. 2: The quantum probability distribution p„{t) for rep- 
resentative simulations is imaged as a function of time (right) . 
The short-time energy-spreading profile is determined by the 
perturbation matrix |W„,noP (left). The number of particles 
is A'' = 50. The upper set is for u = 50, and the lower is for 
u — 5. The strength of the driving is Kd/ Kq = 0.1. For the 
time axis we use dimensionless units r = (-Emax — E(i)t/h, and 
the scaled driving frequency in both simulations is f2 « 0.03. 
The image of the initial level is vertically zoomed, and it has 
the energy e ~ 0.5. The boundaries of the chaotic sea in the 
lower image are indicated by the horizontal dashed lines. 



sical dynamics. But in the chaotic regime one observes 
that the driving is capable of inducing diffusive-like en- 
ergy spreading. This diffusive spreading is restricted to 
the chaotic energy window, and features remarkable cor- 
respondence with the classical simulation. 

Linear response.— The diffusive energy spreading 
in the chaotic regime can be quantified with the time 
evolution of the energy variance. In Fig.[3] we plot the 
time evolution of Ae for both the chaotic and the in- 
tegrable cases. In both we compare the dispersion ob- 
tained under the classical equations of motion for the 
driven system, starting from a micro-canonical ensem- 
ble, to that obtained from quantum evolution from an 
eigenstate with the same energy. As anticipated for dif- 
fusive energy spreading, we observe that in the chaotic 
regime (Ae)^ « 2Dt with diffusion coefficient D (x K^, 
as assumed in the Kubo linear response formula Eq.([5]). 

In Fig|4] we compare the diffusive energy distribution 
that is observed in Fig|2j with the solution of the FPE 
using the diffusion coefficient from Eq.([5|, see [a] 
for technical details. The agreement is good, confirming 
that weakly driven chaotic quantum systems can indeed 
exhibit energy diffusion in regimes realistic for experi- 
mental Bose-Hubbard systems, and solidifying the basis 
of our phenomenological argument for the FPE descrip- 
tion of inter-subsystem equilibration. 

Multiple timescales.— Having established quite 
good quantum-to-classical correspondence in the chaotic 
regime, one wonders whether classical dynamics may in- 
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FIG. 3: The scaled variance as a function of the scaled time 
in the classical (left) and in the quantum (right) simulations. 
The value of the scaled interaction parameter is u = 5 (up- 
per panels) and u — 50 (lower panels). Note the different 
scale of the vertical axis for the quantum vs classical simu- 
lations in panels (c) and (d). Clearly quantum-to-classical 
correspondence fails in the quasi-integrable regime. The val- 
ues oi Kd/Ko are 0.025 (blue), 0.05 (green), 0.075 (red), and 
0.1 (cyan). The driving frequency is as in Fig|2] 
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FIG. 4: (a) The evolving spreading profile p{E), referring 
to the simulation that has been imaged in the lower panel of 
Fig[2] The solid lines (calculation) and the associated symbols 
(simulation) are for r = 149 (narrower), and r — 299 (wider), 
and T = 448 (widest). The lines are based on the numerical 
solution of Eq.([TJ . (b) The density of states g{E) (dotted line) 
and the diffusion coefficient D{E) (solid line) were deduced 
from the diagonalization of the BHH and Eq.([5]l, see for 
technical details. 

dicate features that go beyond simple energy diffusion. 
Fig. [5] illustrates the classical time dependence of e{t), 
and characterizes it by its average value and dispersion. 
Since the phase space has dimension greater than two, 
the possibility of Arnold diffusion guarantees that the 
motion is ergodic within the chaotic sea. This means 
that we can regard different trajectory segments as un- 
correlated pieces of the same infinite time trajectory. If 
we had ergodic motion with a well defined characteristic 
time, all the segments would have the same average and 
dispersion. But this is not what we see: the segments 
have large variation in their dispersion, since they do not 
uniformly fill the whole chaotic region. Rather, the tra- 
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FIG. 5: On the left the scaled energy e{t) as a function 
of time is plotted for a few representative trajectories. The 
right panel displays the average value and the dispersion of e 
within the time interval < t < 40000, for the representative 
trajectories (symbols), as well as for many other trajectories 
(points). Low dispersion values reflect the finite probability 
to encounter sticky motion. The parameters are the same as 
in the lower panel of Figj2j with initial points that have the 
energy e = 0.3. 



jectories contain episodes with long dwell times within 
some sticky regions in phase space, whose existence has 
been confirmed with a Poincare section, see 'a\ If we had 
an unlimited computation power, obviously the expecta- 
tion is to have coincidence of all the points in the right 
panel of Fig[5j But in practice we can address only finite 
time intervals, and therefore the points are scattered over 
a large range. 

Discussion.— We have considered few-mode Bose- 
Hubbard systems as tunably chaotic systems, which in 
chaotic regimes respond generically to weak driving with 
energy diffusion at a rate proportional to the square of 
the driving strength, K^. Consequently we deduce that 
the thermalization of coupled Bosc-Hubbard sub-systems 
can plausibly be described by a phenomenological FPE, 
namely Eq.Q. 

We have also obtained some insight on how detailed 
studies of Bose-Hubbard systems may usefully be com- 
pared to such phenomenological theories. A small sub- 
system can exhibit multiple time scales in its equili- 
bration, because its phase spaces contain sticky regions 
with long dwell times. In principle this may give rise 
to non-Gaussian features, and Levy flight-like deviations 
from strict diffusive behavior. Even if these possibilities 
do not manifest strongly in energy spreading of small, 
driven systems, because the explored phase space vol- 
ume is small, they may perhaps become important in 
multi-component composite systems. It would be impor- 
tant to recognize, then, in interpreting experiments or 
simulations intended to test Eq.Q, that some deviations 
from its diffusive assumptions may not represent errors 
in its depiction of the mesoscopic onset of irreversibility. 



but only the fading traces of microscopic behavior. [a] See supplementary material at [URL] 
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The role of chaos.— The key obstacle for ther- 
malization can be appreciated by considering the com- 
mon paradigm for driven integrable system: the so called 
"kicked rotor" as described by the "standard map" iSlj . 
In the absence of driving the system is integrable. The 
driving amplitude is K. Below a critical value Kc ~ 
0.97 there is no diffusion in energy due to Kolmogorov- 
Arnold-Moser blocking. For somewhat larger values 
D (X {K — Kc)^. Only for strong driving amplitude one 
observes a quasi-linear dependence D cx Kj. 

The dependence of D on the driving amplitude K is 
strikingly different in the case of a driven chaotic system: 
Following the ergodic adiabatic theorem of [M], it has 
been realized [MTO] that a linear response dependence 
D oc shows up for arbitrarily small driving amplitude 
with arbitrarily small driving frequency ("DC limit"). 

In the quantum domain the applicability of linear 
response has been first challenged [TT] and later re- 



analyzed and established |12j using a random matrix 
theory (RMT) and semi-classical perspectives. The ex- 
istence of the underlying classical dynamics is essential 
in order to avoid RMT anomalies that arise beyond the 
regime of 1st order perturbation theory |S2| . 

Strangely enough the semi-classical implied robustness 
of the quantum diffusive behavior has never been verified, 
to the best of our knowledge, for a realistic quantized 
system. More precisely - there are numerous simulations 
in the quantum adiabatic regime where the transitions 
are mainly between neighboring levels. But the regime 
of our interest is different: our interest is in driving in- 
tensities that can be regarded as quantum mechanically 
large, but still semi-classically small. This is the regime 
where the energy landscape can be regarded as a quasi- 
continuum and quantum-to-classical correspondence can 
be expected. 



Obtaining the Einstein relation.— Let us see how Eq.Q is obtained from Eq.([3|. We assume that both systems 
are independently in a canonical state. Accordingly the joint probability distribution is 
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Doing integration by parts on the term that involves d/^D, and noting that dsgi^2 — ±5i,2/3i,2 is a contribution that 
cancels with the (/?i — I32)D term, one observes that we are left with 
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One identifies that this is, up to the inverse temperature factor, merely the canonical average over £), leading to the 
desired result Eq.Q. 



Calculating the diffusion coefficient. The power 
spectrum of W due to the evolution that is generated by 
Hq can be calculated from its matrix elements as follows: 

S{uj) = FourierTransform {W{t)W{Q)) (10) 
= 1^™."!^ 27r(5(a; - (£„-£„)) 

n m 

Here p„ are the occupation probabilities of the T-L^ eigen- 
states. In order to calculate S{uj) for a given micro- 
canonical energy e the practical procedure is to plot the 
smoothed value \W\'^ of the squared elements \Wm,n\'^ as 
a function of e = En along the diagonal {Em — En) = to. 



Then it follows that 

S{uj) = 2Trg{e)\m\ (11) 

If multiplied by the strength of the driving one 
obtains the Fermi-golden-rule expression for the rate of 
transitions due to a monochromatic driving. As implied 
by the Kubo formula Eq. ([5| the diffusion coefficient is 
given by [S3j : 

D{e) ^ (K.nr g{e)\m (12) 

where has implicit dependence on both e and 17 as 
explained above. 



Semiclassical form of the trimer Hamiltonian. 

In a semi-classical context one substitutes ai — y/riie^^^ , 
and defines qi — ipi ~ ipo. Dropping a constant that de- 
pends of the conserved total particle number N, the BHH 
takes the form 

H = -Kf) ^ [noni]^/^ cosgi - [/ [noni-|-noW2+«iW2] 

4 = 1,2 

Expressing no — N—ni—n2 we see that the BHH is the 
quantized version of two coupled degrees of freedoms, 
where the effect of the interaction term (a quadratic func- 
tion of rii and 712), is characterized by the dimensionless 
parameter u — NU/Kq. 

Level statistics and the Brody parameter. 

Given N and K and U, we find the eigen-energies of the 
Hamiltonian Eq. (|6]). In each small energy range we cal- 
culate the mean level spacing, and the distribution P{S) 
of the normalized spacings. Then we fit it to the Brody 
distribution [27j 

Pq{S) = aS"^ exp(-/3S'i+') (13) 

with a = {l + q)(3, and /3 = T^+m2 + q)/l + q)]. Here 
r denotes the Euler gamma function. A Brody param- 
eter value of q = indicates a Poissonian level-spacing 
distribution characteristic of the uncorrelated levels of 
integrable system. By contrast for q = 1 we have the 
Wigner level-spacing distribution, that reflects the level 
repulsion in the case of a quantized chaotic system. Thus, 
by plotting g as a function of e we can map the domain 
of chaotic motion, see Fig[l]3. 

For large u the dynamics is quasi-integrable due to 
self-trapped motion. The dynamics is trivially integrable 
also in the other extreme of very small u. While we 
have also studied this latter region, the results concerning 
the response to driving were similar to the self-trapping 
integrability, and this point is not explicitly discussed in 
order to avoid redundancy. 

Poincare sections.— The left panel of Fig [6] displays 
the Poincare section of a representative classical chaotic 
trajectory. The right panel use the same coordinates for 
plotting a trajectory of the driven system. The sticky 
part of the trajectory is highlighted. 
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FIG. 6: Left panel: The Poincare section of a representative 
classical chaotic trajectory. The parameters are u — 5 and 
e = 0.3. The coordinates qt = (fi — ipo are conjugate to rii, 
with i — 1,2, and the section is at {qi — 52) = 7r/2. Left panel: 
The square labeled trajectory of Fig|5] is illustrated using the 
same coordinates as in the left panel. The points along the 
trajectory that have low energy (e < 0.2) are highlighted by 
dark color: they reside within a sticky region. 
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